Subwavelength modulational instability and plasmon oscillons in nanoparticle arrays 
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We study modulational instability in nonlinear arrays of subwavelength metallic nanoparticles, 
and analyze numerically nonlinear scenarios of the instability development. We demonstrate that 
modulational instability can lead to the formation of regular periodic or quasi-periodic modula- 
tions of the polarization. We reveal that such nonlinear nanoparticle arrays can support long-lived 
standing and moving oscillating nonlinear localized modes - plasmon oscillons. 
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Nonlinearity-induced instabilities are observed in 
many different branches of physics, and they provide 
probably the most dramatic manifestation of strongly 
nonlinear effects that can occur in Nature. Modulational 
instability (MI) in optics manifests itself in a decay of 
broad optical beams (or quasicontinuous wave pulses) 
into optical filaments (or pulse trains) [1], and such ef- 
fects are well documented in both theory and experiment. 
MI is also observed for partially spatially incoherent light 
beams in noninstantaneous nonlinear media with the pat- 
tern formation from noise [2]. It is expected that the 
study of subwavelength nonlinear systems such as metal- 
lic nanowires or nanoparticle arrays may bring many new 
features to the physics of MI and the scenarios of its de- 
velopment, however such effects were never studied be- 
fore. 

Over the past decade, surface plasmon polaritons (or 
plasmons) were suggested as the mean to overcome the 
diffraction limit in optical systems. In particular, by us- 
ing plasmons excited in a chain of resonantly coupled 
metallic nanoparticles [3], one can spatially confine and 
manipulate optical energy over distances much smaller 
than the wavelength. In addition, strong geometric con- 
finement can boost efficiency of nonlinear optical effects, 
including the existence of subwavelength solitons [4]. 

In this Letter, we study modulational instability in 
subwavelength nonlinear systems for an array of opti- 
cally driven metallic nanoparticles [5-8] with a nonlinear 
response. We demonstrate the existence of novel types 
of nonlinear effects in such subwavelength systems never 
discussed before, including the generation of regular or 
quasi-periodic polarization patterns and oscillating local- 
ized modes which can be termed oscillons, in analogy 
with the similar localized modes excited in driven gran- 
ular materials [9] and Newtonian fluids [10]. 

Figure 1 shows the geometry of our problem: a chain of 
identical spherical silver nanoparticles is embedded into a 
fused silica host medium with permittivity su and driven 
by an external optical field with the frequency close to 
the frequency of the surface plasmon resonance of an in- 
dividual particle. We assume that the particle radius 



and distance between the particles are a = 10 nm and 
d = 30 nm, respectively. Ratio a/d satisfies the condi- 
tion a/d < 1/3, so that we can employ the point dipole 
approximation [5]. In the optical spectral range, a linear 
part of silver dielectric constant can be written in a gen- 
eralized Drude form e\ g = — uo^/[uj{ijJ — zz/)], where 
Soo = 4.96, hoj p = 9.54 eV, hv = 0.55 eV [11] (here- 
inafter we accept exp(zcjt) time dependence); whereas 
dispersion of Si02 is neglected since Eh — 2.15 for wave- 
lengths 350-450 nm [12]. Nonlinear dielectric constant 
of silver is = e\ g + x (3) |E^ n) | 2 , where E^ n) is the 
local field inside n-th particle. We keep only cubic sus- 
ceptibility due to spherical symmetry of particles. Cur- 
rently, there is no reliable theoretical models describing 
nonlinear optical response of metal nanoparticles, how- 
ever experimental data shows that x^ depends on many 
factors, including duration and frequency of the external 
excitation as well as particle characteristics themselves 
(metal type and size) [13]. According to the model sug- 
gested in Ref. [14] and confirmed in experiment, 10 nm 
radii Ag spheres possess a remarkably high and purely 
real cubic susceptibility x^ — 3 x 10 -9 esu, in com- 
paring to which the cubic nonlinearity of Si02 is weak 
(~ 10- 15 esu [15]). 

We study nonlinear dynamics of our chain by employ- 
ing the dispersion relation method [16] that allows deriv- 
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FIG. 1. (Color online) Schematic sketch demonstrating ge- 
ometry of the studied problem. Arrows indicate particle po- 
larizations after MI development. 
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ing a system of coupled equations for slowly varying am- 
plitudes of the particle dipole moments. This approach 
is based on the assumption that in the system there are 
small and large time scales, which in our case is fulfilled 
automatically since each particle acts as a resonantly ex- 
cited oscillator with slow (in comparison with the light 
period) inertial response. 

We start with the standard expression for the electric 
dipole moment induced in the n-th particle written for 
Fourier transforms 



a n (u)- 1 p n =-E<£ x) + ^E„, m , 



(i) 



where 



~NL/ 



a n (uj) = eh 



is the electric polarizability of the n-th particle, E^ 63 ^ 
the external electric field acting on n-th particle, 
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is in charge of dipole-dipole interaction between m-th and 
n-th particles, k = Cju/c^/sh, ro is the unit vector point- 
ing from the m-th to the n-th particle. Assuming that 
^( 3 )|En n ^| 2 <C 1 and v/ujq <C 1, we decompose a n (o;) _1 
in the vicinity of the frequency of the surface plasmon 
resonance of an individual particle, ujq = uo p / V^oo + 2^, 
and keep the first-order terms involving time derivatives 
for describing (actually small) broadening of the particle 
polarization spectrum, 
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(2) 



where Auj is the frequency shift from the resonance value. 
Having expressed En^ via p n , we substitute Eq. (2) 
into Eq. (1) and write E n?m in the same order of the 
perturbation theory and obtain the equations, 
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dP n 
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(_ i7 + Q + |p„|2) p- + Gt, m Pm = K, 

- (-i 7 + n + |P„| 2 ) pi + GlmPL = El 

(3) 
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FIG. 2. (Color online) Nonlinear dispersions of (a) longitudi- 
nal and (b) transverse eigenmodes of an infinite chain. Dashed 
curves correspond to the linear limit. Vertical dashed line in 
(b) marks the light line, ko = cjo/cy^. 
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are dimensionless 
slowly varying amplitudes of the particle dipole moments 
and external electric field, respectively, the indices V 
and Y stand for the transverse and longitudinal compo- 
nents with respect to the chain axis, 77 = g 3 ^ g (^) , 
|P„| 2 = \P-\ 2 + \P^\ 2 , 7 = is/(2co ) + (k afe7/(e2+2e h ) 
describes thermal and radiation losses of particles, fco = 
^o/cy^, Q = (uo — ujq) / uoq and r = ujot. Equations (3) 
describe temporal nonlinear dynamics of a chain of metal- 
lic nanoparticles driven by arbitrary external optical field 
with the frequency uo ~ uoq. We stress that the suggested 
model takes into account all particle interactions through 
the dipole fields, and it can be applied to both finite and 
infinite chains, being also extended to higher dimensions. 

First, we consider an infinite chain. For the station- 
ary unbiased linear case, when d/dr = and E^ Al = 0, 
we look for solutions in the form P^ Al ~ exp (—inKd), 
and from Eqs. (3) find well-known dispersion relations 
for transverse and longitudinal eigenmodes of the sys- 
tem [7], shown in Fig. 2. Taking into account nonlinear- 
ity just shifts the dispersion curves along the frequency 
axis. Light line, which for ft « 1 takes the form 
ko = ^o/cy^, divides the eigenmodes into fast (with 
K < ko) and slow (with K > ko) experienced strong and 
weak radiation damping, respectively. Logarithmic sin- 
gularity occurred at K = fco f° r the transverse modes is 
caused by the phase matching between the chain mode 
and the plane wave traveling in the host medium. 

To study MI, we excite the chain by an homogenous 
electric field with one of the two polarizations: (i) E n = 
(Eq , 0) and (ii) E n = (0,£^). In this case, all particle 
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FIG. 3. (Color online) (a) Bifurcation diagram showing bista- 
bility regime and different scenarios of modulation instability 
development for longitudinal excitations, as a function of ft 
and \E\\ 2 . 3 ) Contour map of Ay on the plane (Kd, \Eq\ 2 ) 
at ft = —0.07. Horizontal dashed lines (1) and (2) mark the 
intensities of the external light used in numerical simulations 
of Eq. (3) shown in Figs. 4(a,b), respectively. 



dipole moments remain the same, P^ 11 = Pq'\ and the 
system stationary states can be written as follows 



-i>y + ft - 



£4' 

.7 = 1 



(4) 



where Ai = r] ((&o<f) 2 j -1 — ikodj~ 2 — j -3 ) exp (—ikodj), 
~ 2 + j -3 ) exp (—ikodj). Transition 
has been made via the replace- 
ment \n — m\ = j and taking into account symmetry 
structure of the series. When ft < — Re YlTLi ^t'" — 
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Im^j^zi '" ) 5 the polarization Pq^'" becomes 



a three- valued function of E leading to bistability. 

Next, we analyze linear stability of the stationary 
states with respect to weak spatiotemporal modulations 
and derive the expression for the instability growth rate, 
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where j ±Al = TmZ?=i - 7, = cos(i^j). 

Thus, the initial nonlinear homogenous states (4) become 
unstable provided A^u > 0. The stability depends on the 
external field parameters Eq '" and ft as well as on the 
modulation wavenumber K. 

Next, we consider the case of the longitudinal excita- 
tion in detail. The condition A y = at any K defines 
the boundaries of MI in the plane (ft, |^q| 2 ) snown m 
Fig. 3(a). Interestingly, the middle and upper branches 
in the bistable region of dependency Pq(Eq) also corre- 
spond to MI, but the development inside the bistability 
region cannot be reached because the middle branch is 
unstable, while the system transition from the lower to 
upper branch itself initiates appearance of MI. 

Figure 3(b) shows a contour map of A y in the plane 
(Kd, \E l l \ 2 ) at ft = -0.07. Remarkably, MI takes place 
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FIG. 4. (Color online) Dynamics of Re Pjl obtained by numer- 
ical simulations of Eq. (3) for a finite chain excited longitudi- 
nally with (a) ft = -0.07, |^| 2 = 0.1 x 10" 4 ; (b) ft = -0.07, 
|^| 2 = 4.9 x 10" 4 ; (c) ft = -0.02, |^| 2 = 2.65 x 10" 4 ; and 
(d) ft = -0.09 and |^| 2 = 0.11 x 10" 4 . 



only for slow eigenmodes of the chain. As follows from 
Fig. 3(b), one can manage eigenmode spectrum excited 
during MI growth by varying \Eq\ 2 only. In particular, 
when \Eq | 2 is chosen to be close to the lower or upper 
edge of the MI domain, just one spatial harmonic should 
be excited, with correspondingly Kd c± 0.7 or Kd = it. 

However, the linear stability analysis does not provide 
any information about the subsequent evolution of the 
unstable system, especially when the external field ex- 
cites a broad spectrum of eigenmodes. To analyze those 
scenarios, we perform numerical simulations of Eq. (3) 
for a finite chain (with 100 nanoparticles) at zero initial 
conditions. Edge effects play a role of small perturbations 
needed for generating MI. The amplitude of the homo- 
geneous external field is supposed to be slowly growing 
to the saturation level Eq (which is reached at r « 100) 
lying in the MI zone. 

Characteristic results are summarized in Figs. 4(a-d). 
When E l l crosses the lower edge of the MI domain [de- 
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FIG. 5. (Color online) (a) Bifurcation diagram for transverse 
excitation, (b) Contour map of on the plane of parameters 
(Kd, \Eq\ 2 ) at ft = -0.07. (c) Snapshot of Re at r = 
700 obtained numerically from Eq. (3) with ft = —0.07 and 
|i^j~| 2 = 1.5 x 10 -4 , indicated in (b) by a horizontal dashed 
line. Points joined by dashed lines give a guide for the eye. 



picted by a dashed line (1) in Fig. 3(b)], we observe 
that MI results in the excitation of one eigenmode with 
Kd ~ 0.7 [see Fig. 4(a)], in accord with the prediction 
of the linear stability analysis. The excited eigenmode 
acts as modulation of the initial almost homogenous state 
which becomes unstable. That is why Re P^ tend to be 
predominantly positive. They are just biased by the ex- 
ternal field. 

Figure 4(b) shows the case when the external field of 
larger amplitude excites a wide eigenmode spectrum [in- 
dicated by dashed line (2) in Fig. 3(b)] [17]. Here, MI 
leads to the formation of a stationary higher-order mode 
along with oscillating localized states. Some of them ap- 
peared to be unstable and decay, whereas others remain 
stable. Importantly, such soliton-like localized modes 
may be at rest or they can drift slowly along the chain, 
as shown in Figs. 4 (b,c). We notice that these oscilla- 
tory localized states in a driven chain are very similar to 
spatiotemporal structures termed oscillons observed pre- 
viously in other types of dissipative systems [9, 10], and 
we refer to them as plasmon oscillons. We point out that 
the plasmon oscillons may exist not only in the form of 
solitary states but they also can create patterns, as il- 
lustrated in Fig. 4(d) [17]. We have studied some of the 
properties of such oscillatory states, and the results will 
be published elsewhere. 

Finally, we conduct the similar analysis for the case 



of the transversal excitation. Figure 5(a) shows the cor- 
responding bifurcation diagram. In contrast to the lon- 
gitudinal case, the MI region is fully placed inside the 
bistability domain capturing a part of the lower branch in 
the dependence Pq(Eq). According to the contour map 
of X±_ in the plane (Kd, \Eq | 2 ) shown in Fig. 5(b), the 
spectrum of excited eigenmodes can be tuned by vary- 
ing the value of \Eq | 2 , in analogy with the longitudinal 
case. Nevertheless, the width of spectrum weakly affects 
the scenarios of the MI development. Numerical simula- 
tions of Eq. (3) demonstrate that, independently on the 
value of Eq , the growth of MI results in switching of the 
system from the lower to upper branch in the bistabil- 
ity region of P ± (E L ). In the case of a finite chain, MI 
is accompanied by a pair of switching waves (kinks) at 
the edges which move towards each other, as shown in 
Fig. 5(c) [17]. 

The observation of MI requires high illuminating pow- 
ers, higher than 10 MW/cm 2 , that could cause thermal 
damage to particles. To estimate maximal duration of 
the external laser pulse, we use the results of previous 
studies on the ablation thresholds for gold films [18] pro- 
viding values of 1.6 J/cm 2 and 0.6 J/cm 2 for 1 ns and 
1 ps pulses, respectively. Gold demonstrates stronger 
thermal losses than silver at optical frequencies. That 
is why this data is completely acceptable. Taking into 
account amplification of electric field inside nanoparti- 
cles due to surface plasmon resonance, we come to the 
external threshold intensities of 3.6 MW/cm 2 and 1.3 
GW/cm 2 corresponding to 1 ns and 1 ps pulses, respec- 
tively. Thus, ablation of silver particles will not be crit- 
ical at least till pulse durations of 1 ps. As the charac- 
teristic time of the MI growth is of (X±,\\UJo)~ 1 ~ 10 fs 
that is much less than the maximal pulse duration, all 
predicted effects seem readily observable in experiment. 

In conclusion, we have studied theoretically modu- 
lational instability in arrays of subwavelength metallic 
nanoparticles, and analyzed numerically the development 
of such instabilities beyond the linear approximation. We 
have observed that modulational instability can be en- 
hanced substantially by the geometric confinement, and 
it can lead to the formation of regular periodic or quasi- 
periodic polarization patterns. We have observed the 
generation of long-lived standing and moving oscillating 
nonlinear localized modes in the form of plasmon oscil- 
lons. The experimental observation of the predicted mod- 
ulational instability can provide a prominent approach to 
achieve subwavelength confinement of the optical fields 
guided by plasmonic nanostructures. 
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